package scu.maqiang.possion;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T31;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.mesh.T3Type;
import scu.maqiang.numeric.*;

import java.util.Arrays;

public class TestConvergence2DT31 {
    public static void main(String[] args) {
        ScalarFunc exactFunc = ( xy, label, param) -> Math.cos(Math.PI * xy[0]) * Math.cos(Math.PI * xy[1]);
//        ScalarFunc dxFunc = (xy, label, param) -> -Math.PI * Math.sin(Math.PI * xy[0]) * Math.cos(Math.PI * xy[1]);
//        ScalarFunc dyFunc = (xy, label, param) -> -Math.PI * Math.cos(Math.PI * xy[0]) * Math.sin(Math.PI * xy[1]);
        ScalarFunc source = ( xy, label, param) -> 2*Math.PI*Math.PI*Math.cos(Math.PI * xy[0]) * Math.cos(Math.PI * xy[1]);

        Mesh2T3 rectMesh = new Mesh2T3();
        int n = 8;
        int size = 4;
        double L2Error, L2Error0 = 0, H1Error, H1Error0 = 0;
        double[] L2Conv = new double[n - 1];
        double[] H1Conv = new double[n - 1];
        for(int i = 0; i < n; i++) {
            size *= 2;
            rectMesh.square2D(size, size, T3Type.Cross1);
            FES2T31 fs = new FES2T31(rectMesh);

            int dof = fs.GetNdof();
            SRMatrix A = new SRMatrix(dof);
            double[] RHS = new double[dof];
            double[] x = new double[dof];
            fs.assembleStiff(new double[]{1.0}, BVPType.COMMON, A);
            fs.assembleSource(new ScalarFunc[]{source}, null, BVPType.COMMON, RHS);

            fs.applyBC_MBN(A, Direct.X, 1, 2, 3, 4);
            fs.applyBC_MBN(RHS, Direct.X, new ScalarFunc[]{exactFunc}, null, 1, 2, 3, 4);

            NewIterSSolver solver = new NewIterSSolver(A);
            solver.PCGSSOR(RHS, x, 1.5, 1);

            double[] exactSolution = rectMesh.functionInMesh(exactFunc, null);
            double[] error = MVO.add(exactSolution, -1.0,  x);
            if (i == 0) {
                L2Error0 = fs.computeL2Norm(error);
                H1Error0 = fs.computeH1SemiNorm(error);
            } else {
                L2Error = fs.computeL2Norm(error);
                H1Error = fs.computeH1SemiNorm(error);
                L2Conv[i - 1] = L2Error0 / L2Error / 2;
                H1Conv[i - 1] = H1Error0 / H1Error / 2;
                L2Error0 = L2Error;
                H1Error0 = H1Error;
            }
        }

        System.out.println("L2 Convergence: " + Arrays.toString(L2Conv));
        System.out.println("H1 Convergence: " + Arrays.toString(H1Conv));
    }
}
